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Abstract. We present an analysis method that allows us to estimate 
the Galactic formation of radio pulsar populations based on their ob- 
served properties and our understanding of survey selection effects. More 
importantly, this method allows us to assign a statistical significance to 
such rate estimates and calculate the allowed ranges of values at various 
confidence levels. Here, we apply the method to the question of the double 
neutron star (NS-NS) coalescence rate using the current observed sample, 
and we find calculate the most likely value for the total Galactic coales- 
cence rate to lie in the range 3 — 22 Myr -1 , for different pulsar population 
models. The corresponding range of expected detection rates of NS-NS 
inspiral are (1 — 9) x 10~ 3 yr _1 for the initial LIGO, and 6 — 50 yr _1 for 
the advanced LIGO. Based on this newly developed statistical method, 
we also calculate the probability distribution for the expected number of 
pulsars that could be observed by the Parkes Multibeam survey, when 
acceleration searches will alleviate the effects of Doppler smearing due to 
orbital motions. We suggest that the Parkes survey will probably detect 
1-2 new binary pulsars like PSRs B1913+16 and/or B1534+12. 



1. Introduction 



The detection of the double neutron star (NS-NS) prototype PSR B1913+16 
as a binary pulsar (Hulse & Taylor 1975) and its orbital decay due to emission 
of gravitational waves have inspired a number of quantitative estimates of the 
coalescence rate, 7Z, of NS-NS binaries (Clark et al. 1979; Narayan et al. 1991; 
Phinney 1991; Curran &; Lorimer 1995). Significant interest derives from their 
importance as gravitational-wave sources for the upcoming ground-based laser 
interferometers (such as LIGO). 

We present a newly developed statistical analysis that allows the calculation 
of statistical confidence levels associated with rate estimates. The method can 
be applied to any radio pulsar population. Here, we consider PSR B1913+16 
(Hulse & Taylor 1975) and PSR B1534+12 (Wolszczan 1991). For different 
assumed distributions of pulsar properties (luminosities, Galactic positions), we 
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derive the probability distribution function of the total Galactic coalescence 
rate weighted by the two observed binary systems. The method involves the 
simulation of selection effects inherent in all relevant radio pulsar surveys and 
a Bayesian statistical analysis for the probability distribution of 1Z. The small- 
number bias and the effect of the faint-end of the luminosity function, previously 
identified as the main sources of uncertainty in rate estimates (Kalogera et al. 
2001) are implicitly included in this analysis. We extrapolate the Galactic rate 
to cover the detection volume of LIGO and estimate the most likely detection 
rates of NS-NS inspiral events for the initial and advanced LIGO. Details of this 
work are given in Kim et al. (2003; hereafter KKL). 

In the second part of this paper, we modify our statistical method in a way 
that allows us to calculate the probability distribution for the number of pulsars 
that could be detected by the Parkes Multibeam survey (hereafter PMB survey; 
Lyne et al. 2000; Manchester et al. 2001) PMB, when the effects of Doppler 
smearing due to orbital motions are corrected with acceleration searches. 

2. Pulsar Survey Selection Effects 

For a model pulsar population with a given spatial and luminosity distribution, 
we determine the fraction of the total population which are actually detectable by 
current large-scale pulsar surveys. In order to do this, we calculate the effective 
signal-to-noise ratio for each model pulsar in each survey, and compare this with 
the corresponding detection threshold. Only those pulsars which are nominally 
above the threshold count as being detectable. After performing this process on 
the entire model pulsar population of size AT tot , we are left with a sample of -/V b s 
pulsars that are nominally detectable by the surveys. By repeating this process 
many times, we can determine the probability distribution of iVobs, which we 
then use to constrain the population and coalescence rate of NS-NS binaries. 
More details are given in § 2 of Lorimer et al. (1993) and in KKL. 

Here we discuss in a some detail the the Doppler smearing effect which is sig- 
nificant for recent surveys with relatively long exposures. For binary pulsars, we 
need to take account of the reduction in signal-to-noise ratio due to the Doppler 
shift in period during an observation. For observations of NS-NS binaries, where 
the orbital periods are of the order of 10 hours or less, the apparent pulse period 
can change significantly during a search observation causing the received power 
to be spread over a number of frequency bins in the Fourier domain. As all 
the surveys considered in this analysis search for periodicities in the amplitude 
spectrum of the Fourier transform of the time series, a signal spread over several 
bins can result in a loss of signal-to- noise ratio a. To take account of this effect 
in our survey simulations, we need to multiply the apparent flux density of each 
model pulsar by a "degradation factor", F = ^binary /^control- Significant degra- 
dation occurs, when F < 1. Using an analysis method described in Camilo et 
al. (2000), we calculate the degradation factor for the two pulsars we consider in 
this work. As expected, we find that surveys with the longest integration times 
are most affected by Doppler smearing. For the PMB survey, which has an inte- 
gration time of 35 min, mean values of F are 0.7 and 0.3 for PSR B1913+16 and 
PSR B1534+12 respectively. The greater degradation for PSR B1534+12 is due 
to its mildly eccentric orbit (e ~ 0.3 versus 0.6 for PSR B1913+16) which results 
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in a much more persistent change in apparent pulse period when averaged over 
the entire orbit. In order to improve on the sensitivity to binary pulsars, the 
PMB survey data are now being reprocessed using various algorithms designed 
to account for binary motion during the integration time (Faulkner et al.; these 
proceedings). For the Jodrell Bank and Swinburne surveys (Nicastro et al. 1995; 
Edwards et al. 2001), which both have integration times of order 5 min, we find 
F ~ 0.9 for both systems. For all other surveys, which have significantly shorter 
integration times, no significant degradation is seen, and we take F = 1. 



3. Probability Distribution of Double Neutron Star Coalescence Rates 
3.1. Statistical Method 

As already mentioned, we generate large numbers of "observed" pulsar samples 
by modeling the survey selection effects and applying them to model populations 
of PSR B1913+16-like and PSR B1534+12-like pulsars, separately. For a fixed 
value of iVtot) we use these "observed" samples to calculate their distribution, 
which we find to be very well described by a Poisson distribution: 

\N ohB p -X 

where A =< N -^ s >. With our Monte Carlo simulations we calculate A and find 
it to linearly correlate with N to t (for values in the range 10- 10 4 ): 

A = aJVtot, (2) 

where a is a constant that depends on the properties (space and luminosity 
distributions and pulse period and width) of the Galactic pulsar population. 

For a given A r tot , we calculate the rate using estimates of the associated pul- 
sar beaming correction factor /b and lifetime riif e : TZ = ^r^fh- We adopt values 

discussed in Kalogera et al. (2001): 5.72 and 3.65 x 10 8 yr for PSR B1913+16, 
and 6.45 and 2.9 x 10 9 yr for PSR B1534+12. 

Using Bayes' theorem and the best-fit Poisson distributions, we can cal- 
culate the probability distribution of the total number jV to t of pulsars in the 
Galaxy. Further, using estimates of the associated pulsar beaming correction 
factor /b and lifetime T\if e , we can calculate the distribution function of pulsar 
rates: 1Z = ^^/b- We adopt values discussed in Kalogera et al. (2001): 5.72 

and 3.65 x 10 8 yr for PSR B1913+16, and 6.45 and 2.9x 10 9 yr for PSR B1534+12. 
The probability functions of coalescence rates of pulsars similar to each of the 
observed one are then given by: 

PW = «Ve-^. (3) 
v Jb ' 

We use appropriate variable transformations to then calculate the total rate 
probability distribution: 

(4) 
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where A and B are denned as follows: 




1913 



and B = 



ariife 



fb 



) 



1534 



(5) 



Having calculated the probability distribution of the Galactic coalescence rate, 
we can take one step further and also calculate ranges of values for the rate T^tot 
at various confidence levels (CL). The lower (7£ a ) and upper (TZh) limits to these 
ranges are determined by the following conditions: 



Finally, we can calculate the detection rate for LIGO, TZ^et, defined by 



where e is the scaling factor (based on the blue luminosity density of the nearby 
universe) derived to be ~ 10~ 2 Mpc~ 3 (for details see Kalogera et al. 2001). 
Vdct is the detection volume defined as a sphere with a radius equals to the 
maximum detection distance D max for the initial (~ 20Mpc) and advanced 
LIGO (~ 350 Mpc; Finn 2001). 

3.2. Results 

We have chosen one of our pulsar population models to be our reference model 
based on the results presented by Cordes & Chernoff (1997). For this model, we 
find the most likely value of N tot , to be ~ 390 pulsars for the "PSR B1913+16- 
like" population, and ~ 350 pulsars for the "PSR B1534+12-like" population. 
Using eq. (4-5), we evaluate the total Galactic coalescence rate of NS-NS bi- 
naries for this reference case. The most likely value of the coalescence rate 
is IZpeak— 8Myr _1 and the ranges at different statistical confidence levels are: 
~ 3 - 20 Myr- 1 at 68%, ~ 1 - 30 Myr" 1 at 95%, and ~ 0.7 - 40 Myr" 1 at 99%. 
Also, the most likely values of detection rates, which correspond to 7£p Ga k are 
~ 3 x 10~ 3 yr _1 and ~ 18 yr -1 , for the initial and advanced LIGO. 

In Fig. 1, P(Htot) along with P(TZiqi^) and P(TZi^) are shown for the 
reference model. It is evident that the total rate distribution is dominated by 
that of PSR B1913+16. This is due to the fact that we calculate the two rate 
contributions having relaxed the constraint that pulsars have luminosities equal 
to that of the observed pulsar, and instead allowing for the full range in lumi- 
nosity. In this case any differences in the two separate rate contributions depend 
only on differences in pulse periods, and widths. Given that the latter are rather 
small, it makes sense that, for example, the most likely values of N tot for the two 
pulsars come out to be very similar (e.g. ~ 390 and ~ 350, for PSR B1913+16 
and PSR B1534+12, respectively, in the reference model). Consequently any 
difference in the rate contributions from the two populations is due to the differ- 
ence in lifetimes (about a factor of 10) for the two observed pulsars (note that 
the two do not only have similar N to t estimates, but also similar beaming cor- 
rection factors). Since the lifetime estimate for PSR B1913+16 is much smaller, 
the total rate distribution is dominated by its contribution. 




P(H tot )dH tot = CL and P(K a ) = P{U h ). 



(6) 



T^det — e ^-totKjet, 



(7) 



APS Conf. Ser. Style 



5 



22 
-t-> 

13 
to 

-O 

o 
t- 
a. 



0.01 0.1 1 10 100 

Coalescence Rate (Myr -1 ) 

Figure 1. The probability distribution function of coalescence rates 
in both a logarithmic and a linear scale (small panel). The solid line 
represents PiJZtot) and the long and short dashed lines represent P(TZ) 
for PSR B1913+16-like and PSR B1534+12-like populations, respec- 
tively. The dotted lines indicate the confidence levels for P(7^ Q t)- 

We have found that there are strong correlations between the peak value 
of the total Galactic coalescence rate, 7^ pe ak, and the cut-off luminosity, L m i n , 
and its power index, p. As seen in the Fig. 2, 7^ pC ak increases rapidly with 
decreasing L m i n (left panel) or with increasing p (right panel). The peak rate 
does not show any strong dependence on scale lengths of the spatial distribution 
(either i?o or Zq), except for rather extreme cases (e.g. i?o < 3 kpc). Hence, the 
most important model parameter seems to be the slope and low-end cut-off of 
the luminosity function. We find that the most likely values for the rates are 
in the range ~ 3-22 Myr -1 for all models with luminosity-function parameters 
consistent with pulsar observations at 68% confidence level (Cordes & Chernoff 
1997). 



4. Probability Distribution of ./V obs for Parkes Multibeam Survey 
4.1. Statistical Method 

As mentioned in §2.1, because of the long integration time, the signal-to-noise 
ratio for the PMB survey is severely reduced by Doppler smearing due to the 
pulsars' orbital motion. Acceleration searches in the reanalysis by Faulkner 
et al. (these proceedings) promise to alleviate this reduction in the near future. 
Here we calculate the probability distribution of the number of pulsars -/V bs that 
could be detected with the PMB survey, assuming that the reduction in flux due 
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Figure 2. Left panel: The correlation between 7£ pca k and the cut-off 
luminosity L m i n for different power indices p of the luminosity distri- 
bution function. Right panel: The correlation between 7£ pe ak and the 
power index of the luminosity distribution function p. 



to Doppler smearing is corrected perfectly. Using this distribution P(iV b s ), we 
calculate the average value of N ^ s , < N Q ^ S >pmb> for the PMB survey. 

We first calculate N Q ^ S (detected by PMB survey) for a range of total num- 
ber of pulsars in our model galaxy and determine the slope, opmb using eq. 
(2)(§3.1). We then use P(N tot ) (in practice -P(A)) derived in §3, to calculate 
P{N bs), for each type of pulsar: 

P{N ohs ) = J J P(iVobs;ApMB)P(ApMB) dApMB, (8) 

where P(N ohs ; A PM b) is given in eq. (1). Defining (3 = P(A PM b) can be 

derived from P (N tot ): 

^(ApMB)=/? 2 ApMBe- /3ApMB . (9) 

Since a's are different for each type of pulsar, f3 is determined for each pulsar, 
separately. The normalized probability distribution of iV bs for each type of 
pulsar, P(A/obs), is then calculated as: 

P(N , ) - ^ + 1) (10) 

[ ° bsj " (l+/?)2(l +/ J)iW ^ 

Once we have P(A r b s )i9i3 and P(N b s )i534, we can calculate the combined 
P(N ohS;1913 + iVobs,i534)- We define N + = iV obSi i9i3 + iVobs,i534 and iV_ = 
^obs,i9i3 — ^obs, 15345 where all variables are integers. The joint probability distri- 
bution function for observing either PSR B1913+16-like or PSR B1534+12-like 
pulsar, P(iV + ,iV_) is given by: 

P(N AM r ^ I 2 (N + + N. + 2)(N + -N. + 2) 

1 + ' ^ b(l +/3i)(l +(3 2 )\ (l+0 1 )m+*-)/2(l + l3 2 ){N+-N-)/2- 1 ' 



APS Conf. Ser. Style 7 

0.4 I 1 1 1 1 1 1 1 1 1 



0.3 - 



2°- 2 - 



0.1 - 



o I ' ' ' ' 1 1 1 1 1 

O 1 23456789 10 
N + 

Figure 3. The probability distribution function P(-/V bs) f° r the PMB 
survey in the absence of any signal-to-noise ratio due to orbital mo- 
tions. Results are shown for our reference model, but they are highly 
insensitive to the model parameters. 

Finally, we calculate P(N+) by the summation of P(N + ,N-) over N-. 

P(N + ) = Y / P(N + ,N_), (12) 

N_ 

where iV_ lies in the range [— N + ,N + ] with an increment of 2. 
4.2. Results 

In Fig. 3, we show the distribution function P(N + ) that could be detected by the 
PMB survey for our reference pulsar population model. We have calculated the 
distributions and average predicted values, for other models with different cut-off 
luminosity L m \ n and the power index of the luminosity function, p (KKL). For 
7 different models with luminosity-function parameters within a 68% confidence 
level (Cordes Sz Chernoff 1997), we find values of < N+ > in the range 1.35-1.5 
(1.4 for our reference model). We conclude that, if orbital motion does not affect 
the signal-to-noise ratio, then the PMB survey could be expected to detect 1-2 
new binary pulsars with pulse profile and orbital properties similar to either 
PSR B1913+16 or PSR B1534+12 in the PMB survey. 

5. Discussion 

We have recently developed a new method for estimating the total number of 
pulsars in our Galaxy and have applied it to the calculation of the coalescence 
rate of double neutron star systems in the Galactic field (for more details see 
KKL). Here, we extend this method to obtain a prediction for the average num- 
ber of observed pulsars that the PMB survey could detect when acceleration 
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searches are used to correct for the Doppler smearing due to orbital motions. 
The modeling of pulsar survey selection effects is formulated in a "forward" 
way, by populating the Galaxy with model pulsar populations and calculating 
the likelihood of the real observed sample. This is in contrast to the "inverse" 
way of the calculation of scale factors used in previous studies. 

We note that this method could be further extended to account for distri- 
butions of pulsar populations in pulse periods, widths, and orbital periods. It is 
important to note that both our rate estimates and the predictions for detections 
from the PMB survey do not apply to binary pulsars that are significantly dif- 
ferent from with such properties that are significantly different PSRs B1913+16 
and B1534+12 in terms of pulse shapes and orbital properties. 

Most importantly the method can be applied to any type of pulsar popula- 
tion with appropriate modifications of the modeling of survey selection effects. 
Currently we are working on assessing the contribution of double neutron stars 
formed in globular clusters as well as the formation rate of binary pulsars with 
white dwarf companions that are important for gravitational-wave detection by 
LISA, the space-based interferometer planned by NASA and ESA for the end of 
this decade. 
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